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SUMMARY 

A computer code is under development whereby the thin-layer Reynold3-averaged Navier-Stokes 
equations are to be applied to realistic fighter-aircraft configurations. This transonic 
Navier-Stokes code (TNS) utilizes a zonal approach in order to treat complex geometries and satisfy 
in-core computer memory constraints. The zonal approach has been applied to isolated wing geometries 
in order to facilitate code development. Part 1 of this paper addresses the TNS finite-difference 
algorithm, zonal methodology, and code validation with experimental data. Part 2 of this paper 
addresses some numerical issues such as code robustness, efficiency, and accuracy at high angles of 
attack. Special f ree-3tream-preserving metrics proved an effective way to treat H-mesh singularities 
over a large range of severe flow conditions, including strong leading-edge flow gradients, massive 
shock-induced separation, and 3tall. Futhermore, lift and drag coefficients have been computed for a 
wing up through . Numerical oil flo w patterns and particle trajectories are presented both for 

subcritical and transonic flow. These flow simulations are rich with complex separated flow physics 
and demonstrate the efficiency and robustness of the zonal approach. 

1. Introduction 

Computational fluid dynamics (CFD) has rapidly evolved over the past few years from simulations 
involving simple two-dimensional invi3cid flows about airfoils to three-dimensional wing and 
wing/body computations using the more complete Navier-Stokes equations. This has been attributed to. 
significant improvements in computer speed and memory capabilities, as well as similar advances in 
the development of efficient numerical algorithms 1,2 . Consequently, it is now feasible to undertake 
the task of simulating transonic viscous flow about realistic aircraft configurations with the 
Navier-Stokes equations. A computer code is currently under development at NASA Ames Research Center 
whereby the thin-layer Reynolds-averaged Navier-Stokes equations are being applied to a complete 
fighter-aircraft configuration. Thi3 transonic Navier-Stokes code (TNS) utilizes the Numerical 
Aerodynamic Simulator (NAS) supercomputers at Ames such as the Cray X-MP and Cray 2. The TNS code 
will eventually be used to simulate flight conditions spanning the entire flight envelope of the 
aircraft, including flight conditions at . 

The generation of a single computational grid about a realistic fighter configuration with 
suitable grid clustering is a formidable task. Even if this were accomplished, the grid would 
exceed the current in-core memory limitations of the Cray X-MP computer. In order to overcome this 
difficulty, a zonal approach has been adopted whereby the flow field is subdivided into smaller 
zones. This simplifies the grid generation procedure because existing grid generation methods can 
be applied to these smaller, simpler zones. Furthermore, the zonal approach also combats the memory 
constraint of the Cray X-MP because only one zone and it3 temporal solution need reside in-core at 
any one time. The remaining zones and their temporal solutions reside on the Cray SSD, an efficient 
mass storage device. Thus a zone is brought into core when needed and transferred to the SSD when it 
is not. This transfer process only takes about l| percent of the total CPU time, 

A natural first step in the development of the TNS code is to apply the zonal approach to an 
isolated wing. This permits development and validation of the zonal approach on the simpler wing 
geometry. Moreover, a rich variety of viscous flow physics can also be studied, including highly 
three-dimensional vortical and separated flow. Transonic Navier-Stokes wing solutions are presented 
up through Ci m:lx using the TNS code. The motivation for computing these flows was to determine 
if the TNS code was robust enough to simulate high-angle-of -attack massively separated flow with 
the current zonal topology, and provide a first attempt at computing C[, mar with the Navier-Stokes 
equations . 

The TNS code consists of three major components. The first is the grid generation program. 

There is an option to use the elliptic method of Sorenson 3 or the parabolic method of Edwards 4 . The 
computations presented in this paper use the parabolic grid generator. Once a coarse outer grid is 
generated for the wing, a grid zoning program subdivides the flow field into four zones and refines 
the mesh near the wing surface. The grid generation and zoning programs are coded in a general 
manner and can accomodate a large variety of wing geometries and mesh clustering. Finally, the flow 
solver is applied to the zonal topology and a numerical solution i3 obtained. Present calculations 
utilize four zones totaling 150,000 grid points. 

This is the second part of a two-part paper. Part 1, authored by Flores et al , provides a 
detailed description of the TNS code, including the zonal approach and code validation in comparison 
with experimental data. Part 2 is primarily intended to address numerical issues such as code 
robustness, efficiency, and accuracy, including an effective treatment of coordinate singularities. 
The governing conservation laws and the finite-difference equations (FDE) U3ed to numerically 
integrate these equations are presented in Section 2. For the sake of continuity between Parts 1 and 
2. a brief description of the zonal approach is provided in Section 3. In Section 4, a description 
of special numerical metrics that effectively treat H-mesh coordinate singularities is presented. 
Computational results, obtained for the isolated wing at angles of attack ranging from 0° to 15° are 
presented in Section 5, while concluding remarks are made in Section 6. 

2. GOVERNING EQUATIONS 

The thin-layer Reynolds-averaged Navier-Stokes equations are solved in strong conservation-law 
form and written in generalized curvilinear coordinates as 

d r Q + d ( E + d„F + d,G = Re~ l d^S 


( 1 ) 
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expression the sealed contravariant velocity components are defined by 


In the above 


U = it + €*t» + iv v + ( 2a ) 

V — rft + tj x u + riyV + rj*u/ ( 26 ) 

W = ft + f*« + + ; z w (2c) 


and the pressure is related to the other flow variables by the perfect gas law 


P = (7 ~ 1) 



\p{u 2 + v 2 + w 2 ) 



The transformation metrics are defined as 


Zz = - Zr,y<), Vz = - Vf 2 ?)- (z - J{y^r, - ft = “^rf* - yri v - *r f* 

iy = J(z„x ( - x,« ? ), riy = J(x^2 f - Z£Z f ), ? y = ./fox,, ~ x^z,), 7t = -Xr»>* - yrty - Zrriz 

f* = J(x^ - y„X ; ), h* = -/fox f “ i«y ? ), f* = -/foy,, - y«^), ft = -X T (x - PrCy - Zrf* (4) 

and the Jacobian by 

J~ l = zj(y„z ( ~ z vSAf) + z-ifoy? ~ PfZ?) + x t foz„ - z ( y n ) (5) 

A choice of two implicit approximate-factorization algorithms are available in TNS to integrate 
Eq. (1). The first option is the Beam-Warming algorithm 3 given by 


(/ + h6 ( A n - hDi\ ( )(I + k6 n B n - fr£> t |„)(/ + M f C n - hR e - l S i J- l M n J - hDi\ s )AQ n = R n (6) 

where ^ ^ ^ ^ „ 

R n = -h{S i E n + 6 n F n + S i G n - Re- l ^S n + D t Q n ) 

In the above expression A. B. C, and M, are the Jacobian matrices of E. F. G. and 5. 
respectively. The above method is second-order accurate in time (trapezoidal rule) when h = |A t and 
fir3t-order accurate (Euler implicit) when h — At. The spatial operators use central differencing 
throughout, so fourth-order explicit (Z? e ) and second-order implicit (£>,•) numerical dissipation 
terms are added in order to damp any high frequency errors. This method requires the solution 
of a block tridiagonal system of equations which is computationally costly. For more details of 
this algorithm see Pulliam and Steger 8 . The second option is a diagonal form of the Beam-Warming 
algorithm developed by Pulliam and Chaus3e 7,8 . It has the form 

T ( {I + h6( A? + hDi\ ( )N{I + hS„ A" + hDi\ v )P{I + kS { A? + hD,\ t )T^ l AQ n = R n (7) 


'where is a diagonal matrix consisting of the eigenvalues of A, and so on. The diagonal algorithm. 
Eq. (7). uses fourth-order implicit dissipation, is first-order accurate in time, and only requires 
the solution of scalar pentadiagonal equations. The right hand side of Eq. (7) is the same as 
the right hand side of the Beam-Warming algorithm, Eq. (6). The diagonal algorithm is used for 
steady-state computations or first-order time-accurate integration. 

Convergence to a steady-state solution can be accelerated by using the spatially varying time 
step of Srinivasan et al 9 given by 


h = 


Atp 

1 + \fj 


( 8 ) 


where At 0 is a user specified constant. Flores 10 , using the diagonal algorithm with variable time 
step, reported convergence rates 40 times those obtained using the Beam-Warming algorithm with 
constant time step. 


3. ZONAL APPROACH 

The computational regime about an isolated wing is subdivided into four zones, or grids, with a 
symmetry plane at the wing root. A symmetry-plane view of the four zones is shown in Fig. (1). The 
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coarse outer grid, zone 1, is supplied from the grid-generation program and has a cut-out region that 
is occupied by three finer grids supplied from the grid-zoning program. Throughout the rest of thi3 
paper terms such as zone, grid and block are used interchangeably. 

The grid-zoning program takes block 1 and fills the empty region with block 2, [Fig. (1)]. Like 
block 1, block 2 also has a cut-out region about the wing. The viscous terms are turned off in 
blocks 1 and 2 because they are considered negligible in these zones . Finally blocks 3 and 4 are 
formed in the zoning program with a high degree of clustering in the direction normal to the wing 
surface to support the viscous terms in the boundary layer. Block 3 is on top of the wing and block 
4 is below the wing. 

All four zones overlap their adjacent zones by one or two cell widths, with the exception of the 
interface between blocks 3 and 4. The solution is advanced to the next iteration level (or time 
level for time-accurate computations) by bringing that particular zone into core memory from the Cray 
SSD. Boundary data is supplied from the overlap of adjacent zones. Once the flow has been updated 
in a zone it is transferred to the SSD and the next zone is brought into core memory. Thi3 process 
continues until all zones have been updated. The transfer of boundary data between blocks 2 and 
3 is direct injection and therefore fully conservative. The same is true between blocks 2 and 4. 

The solution at the interface between blocks 3 and 4 is updated by using simple averages. For more 
details of the zonal approach and the TNS code 3ee Holst et al u , Flores 10 , or Part 1 of this paper. 

4. CONSISTENT METRICS 

4.1 Motivation 

The zonal grid topology described in the previous section utilizes Cartesian-like grids to 
simplify the zonal interfacing and maintain flow conservation at shocks. One drawback to this 
approach, however, is the resulting H-me3h singularity at the wing leading edge. This occurs at the 
interface between blocks 3 and 4, as shown Fig. (2). It is important to be able to properly treat 
severe coordinate singularities because they naturally arise in realistic aircraft configurations. 

Some of the early computations with the TNS code were performed at low to moderate angles of 
attack and used simple central differencing of the metrics. These relatively simple test cases 
appeared to have no difficulty at the wing leading edge. In order to control metric truncation 
errors and insure uniform flow as an exact solution of the finite-difference equations, the 
free-stream residual wa3 subtracted from the right side of Eq. (7). However, for M <*, = 0.5 and 
a = 10°, large oscillations in the flow variables developed at the H-mesh singularity as indicated by 
the wing root section Mach number contours given in Fig. (3). It was not possible to obtain a fully 
converged solution. A nominal two order magnitude drop in the L 2 norm of the residual wa3 possible 
only after adding large amounts of numerical dissipation near the leading edge. This was clearly 
unacceptable and would prevent high-angle-of -attack simulations. This problem was fixed by using 
free-stream preserving metrics as described by Pulliam and Steger 1 2 3 . 

In their paper, Pulliam and Steger used special numerical metrics that insured uniform flow was 
an exact solution of the FDE or the free-stream subtraction described above. In the cases treated 
in lief. 6, Pulliam and Steger found no appreciable difference between the two solutions or in their 
convergence rates, and therefore preferred the latter approach. However, tneir grids were relatively 
smooth and did not possess severe coordinate singularities 3uch as the H-type. Once the free-stream 
preserving metrics were implemented in the TNS code, there was no difficulty in rapidly converging 
to a steady state, and no additional dissipation was necessary at the leading edge. The resulting, 
improved, Mach number contours are shown in Figs. (4-5). There is a nominal amount of distortion at 
the interface of blocks 3 and 4 because the flow variables are obtained there using simple averages 
rather than from the governing equations. The norm of the residuals for all four blocks are shown 
in Fig. (6). A three-order drop in the residuals of all four blocks was obtained in 700 iterations. 
This criterion for convergence is usually sufficient for plottable accuracy. 

4.2 Theory 

The use of free-stream preserving metrics in the present finite-difference algorithm has been 
shown to be effective in treating the H-mesh singularity. The necessary conditions the numerical 
metrics must satisfy in order for uniform flow to be an exact solution of the FDE is now presented. 

Upon imposing uniform flow on Eq. (1) and assuming a stationary grid results in 

d(Eoo + drjFoo + dfGca = 0 (9) 

Notice that there is no viscous term in Eq. (9) since S 0 „ = 0. This equation reduces to the following 
three metric relations 

(*),♦($).♦(*),- 

It i3 easily seen with the aid of Eq. (4) that the above relations are exactly satisfied because 
derivative operators obey the two conditions: 

1) Commutative law, i.e. d^d n = d n dc 

2) Chain rule, e.g. dc[uv) = udcv + vdeu. 

Mo3t finite-difference operators satisfy condition one (e.g. central differences); however 
condition two is difficult to satisfy numerically. It is possible, however, to evaluate the metrics 
in a manner such that truncation errors combine to exactly satisfy Eqs . (10). The precise manner in 


(10a) 

(106) 

(10c) 
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which to evaluate the numerical metrics depends on the finite-difference algorithm. For the present 
central-difference algorithm, a consistent set of metrics satisfying Eqs. (10) are given by 

~J = {n t 8„y)(n n 6 t z) - ~ = ((I'Sfzlit i&y) - (/f t %)(Mf<$ f *), j = (iz^y^iz^z) - (iz n Scz)(ixiS,,y) 

y = (it t S„z)(n v S ( x) - (n'SnX^tivSfz), y = (/i { 6 ( x)(ft ( 6 { z) - (^ f S ( z)[n ( S ( *). y = (^8 ( z){n ( 8^x) - [n v S ( x){ix ( S v z) 

y = - (n t 6„y){n,6 t x), -j = {^8 ( y){n ( S^x) - (M«*)(/*e<f»)> y = (M«*)(/*<M) - {y. r ,& i y)(fi e J r ,x) 

(11) 

where ii^,n v , and are simple averaging operators in the £.y, and f directions, respectively. 
Evaluating metrics at boundaries with one-sided differences does not alter the free-stream preserving 
property. 

Equations (10), which are exact relations for the transformed differential equations, have a 
simple geometrical interpretation, quantities such as Z x /J. f) z /J and $ X /J represent the x component 
of area of constant £. constant r) and constant £ surfaces respectively. Thus the divergence-free 
relations, Eqs. (10), represent the necessary condition that the vector surface area of a closed 
volume must sum to zero. The relationship between physical areas and computational areas is given by 

{ dy dz\ / rj x ( x \ /dr) d{\ 

dA = \ dx dz = - f y T) y d? (12) 

\dxdy ) J \(.z f) x ) \di dr) ) 

One can also obtain Eqs. (10) by integrating the above expression about a closed surface and applying 
the divergence theorem. Equations (10) also provide a connection between finite-difference and 
finite -volume methods. 

5. RESULTS AND DISCUSSION 

The numerical wing solutions presented in this section are the result of a performance evaluation 
of the TNS code at high angles of attack. There was 3ome concern regarding code robustness and 
stability because of the H-mesh leading edge singularity and the zonal boundary updating procedure 
(see Figs. 1-2). These solutions are intended to demonstrate that the TNS code, with the consistent 
metrics (Section 4), has the same robustness, efficiency, and stability characteristics at high 
angles of attack as has been observed at low to moderate angles of attack (see Part 1). These flow 
simulations include large leading-edge flow gradients, massive separation, and computations at Ct mul . 

A generic low aspect ratio wing consisting of a NACA 0012 airfoil cross section was selected for 
this investigation. A perspective view- of the wing is shown in Pig. (7). This untapered, untwisted, 
wing has a 20° leading-edge .sweep angle, and an aspect ratio of three. Significant three-dimensional 
effects are encountered due to the low aspect ratio. A symmetry-plane boundary condition is imposed 
at the the wing root section to reduce computational time and computer memory requirements. There 
are a total of 150,000 grid points in all four zones. The grid spacing normal to the wing surface 
varies between 10 -s 10~ 5 chords. 

The diagonal algorithm with spatially-varying time step, Eqs. (7-8), is used to obtain 

steady-state solutions. The Baldwin- Lomax 12 algebraic turbulence model was also used because it is 

efficient and adequate for this present numerical study. Turbulent transition i3 imposed at the wing 
leading edge, and the free-stream Reynolds number is 8 million based on the wing root chord. 

The first set of solutions were obtained at a free stream Mach number Moq = 0.5. A comparison 

between TNS lift and drag coefficients and those obtained with a method based on the full potential 
equations 13 is shown in Figs. (8-9). These flow simulations indicate that the flow is subcritical 
due to the effects of three-dimensional relief and wing sweep, Sonic flow is achieved only at 
maximum lift (a = 13.5°). The lift coefficient obtained with the TNS code i3 in good agreement 
with the full potential result in the low-angle-of -attack range but differs significantly in the 
high-angle-of -attack regime. The TNS code predicts maximum lift at a = 13.5° while the full potential 
lift coefficient continues in a linear fashion. The drag coefficient exhibits the usual quadratic 
variation with angle of attack (Co ~ C 2 ~ a 2 ) together with the large drag rise at stall. The 
convergence history of the above solutions is similar to the convergence history shown in Fig. (6). 

A converged solution usually requires 700 iterations or 55 minutes of Cray X-MP time. 

The NACA 0012 wing geometry described above was used in Part 1 of this paper to compare 
compressible TNS flow simulations with available low-angle-of -attack experimental data. The 
computations reported in thi3 second part of the paper are an extension of the low-angle-of -attack 
computations of Part 1 to the high-angle-of -attack regime. However, compressible experimental 
data at large angles of attack was not available for comparison with computation. Unfortunately, 
there seems to be a lack of good force, moment and surface pressure data up through maximum lift in 
the literature. Although, force and moment data provide a way to asse33 a codes global accuracy, 
extensive surface pressure data is necessary for validating the details of the flow simulation. 

Particle trajectories for a = 15° are shown in Figs. (10-11). Figure (10) is a perspective view 
from above the wing and looking downstream toward the wing leading edge. The particles are released 
along the wing leading edge and wing tip. This massively separated steady flow exhibits a-induced 
separation, i.e. separation due solely to angle of attack, as apposed to shock-induced separation. 

The vortical structure of the separated region is evident. The particle trajectories emanating from 
the wing tip indicate a wing-tip vortex. Figure (11) is an end view of the wing looking inboard from 
the wing tip. The region of separation extends across zonal boundaries in a smooth manner. This 
solution seem3 to be on the verge of going unsteady and required twice as many- iterations to achieve 
convergence. The spatially varying time step tends to inhibit unsteadyness . When a constant time 
step was used, a cyclic variation' in the residuals was noted. 

All of these solutions where obtained with relatively little difficulty, provided consistent 
metrics were used. Without these special free-stream-preserving metrics, solutions at angle of 
attack greater than 5° would not converge. These high-angle-of -attack solutions also demonstrate the 
robustness of the zonal approach. 
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Transonic wing solutions were also obtained for Mco = 0.8. The TNS Ci vs a variation is 
compared with the previous subcritical case in Fig. (12). Maximum lift occurs at a much lower angle 
of attack in this transonic flow, a = 6° , due to shock-induced separation. The convergence histories 
of the residuals and lift coefficients are shown in Figs. (13-14). These transonic solutions were 
relatively easy to obtain, requiring 50 minutes per solution. 

The wing-tip vortex for the maximum lift condition, a = 6°, is shown in Fig. (15). A very 
interesting simulated oil-flow pattern on the upper surface of the wing i3 shown in Fig. (16). 

These are particle trajectories constrained to the next coordinate surface above the wing (due to 
the no-slip condition on the wing) . Notice the saddle and nodal point singularities on the upper 
surface of the wing. These critical points form a stable topological configuration as defined in 
Ref. (14). There is a major separation line extending over most of the wing span and is followed by 
a reattachment line (shown by a dashed line) a short distance downstream. A second 3mall separation 
region is also evident near the trailing edge in the vicinity of the wing tip. 

The extent of separation in these 3Ubcrltical and transonic cases tend to be underpredicted. 

This became evident in Part 1 where comparisons were made between simulated and experimental surface 
oil-flow patterns. There are several reasons for this discrepancy. First, better grid resolution 
is probably required to improve the accuracy of the computation. Grid refinement studies using a 
million point wing grid on the Ames Cray 2 are currently in progress. Second, an improved numerical 
dissipation model is required. Blended fourth-order and 3econd-order smoothing i3 used in Eq. (7) 
as described by Pulliam 15 . This blended smoothing has very good shock-capturing characteristics 
on a grid suitable for the Euler equations. However, the dissipation coefficient varies as C>(xy) 
and, on fine-3paced grids needed for viscous computations, the numerical dissipation can be as large 
as the physical dissipation. Finally, an improved and efficient turbulence model i3 required that 
can adequately model three-dimensional shock-induced separation. The first two points are easily 
achievable in the near future, the latter still remains uncertain. 

6. CONCLUDING REMARKS 

Three-dimensional transonic Navier-Stokes wing computations have been presented within the 
framework of a zonal approach. Consistent frrae-stream preserving metrics were found to be effective 
in treating the leading-edge H-mesh singularity for a variety of subcritical and transonic flow 
conditions, including maximum lift. Lift and drag-coefficient variation with angle of attack were 
presented up through , together with particle trajectories and simulated oil flow patterns 

that identify extensive flow separation and critical points. The TNS code is capable of simulating 
high-angle-of -attack aerodynamics with the zonal approach in an efficient manner, requiring about 
45-60 minutes of Cray X-MP time per solution on a 150,000 point mesh. Three areas are suggested 
for improving the quantitative capabilities of the TNS code, i.e., 1) Grid refinement, 2) improved 
numerical dissipation model that reduces the smoothing within the boundary layer and, 3) an improved, 
efficient turbulence model that treats shock-induced separation. Good experimental force, moment and 
detailed surface pressure data are still needed up through Ct ma , in order to validate compressible 
CFD codes in the high-angle-of -attack range. The robustness and efficiency of the TNS code makes 
it a good candidate as a base code for developing and evaluating new turbulence models capable of 
treating shock-induced separation. The TNS code can al30 be used as a design tool to improve wing 
efficiency. 
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Fig. (1) Symmetry plane view of four zone grid 
topology. 



Fig. (2) H-mesh singularity at wing leading edge, 
(blocks 3 and 4). 




Fig. (3) Symmetry plane Mach number contours 
with central difference metrics; = 0.5, a = 

10.0 °, Re = 8 x 10®. 


Fig. (4) Symmetry plane Mach number contours 
with new metrics; = 0.5, a = 10.0 °,Re = 

8 x 10 6 . 



Fig. (5) Symmetry plane Mach number contours 
with new metrics; M a 0 = 0.5, a = 10.0 °,Re = 
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Fig. (6) L? norm convergence history of residuals for 
all four blocks; M~, = 0.5, a = 10.0°,fj€ = 8 x 10". 












Fig- (7) Perspective view of NACA 


0012 wing. 
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Fig. (10) Perspective view nf „ *■ , 

0Ver a stalled wing- M _ / Particle trajectories 
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Fig. (12) TNS Lift coefficient comparison between 
subcritical (M x — 0.5) and transonic (M x = 0.8) 
cases. 


o 



Fig. (13) L 2 norm convergence history of residuals 
for all four blocks; = 0.8, a = 6.0°, Re — 8 x 
10 s . 


o 


in 



Iteration 



Fig. (14) Lift coefficient convergence history; 
0.8, a = 6.0°, Re = 8 x 10 s . 


Fig. (15) Particle trajectories of the wing tip vortex 
at maximum lift; M x = 0.8, a = 6.0°, Re = 8 x 10 s . 
The wing planform near the wing tip is shaded. 
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Fig. (16) Numerical oil flow pattern on upper wing 
surface; = 0.8, a = 6.0°, Re = 8 x 10 s . 

(N — nodal critical point, S — saddle critical point). 
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